#include "atmoslayer.h"
#include "mie_c_src/mie.h"
#include <QFile>
#include <QByteArray>
//#include <QDataStream>
#include <QTextStream>
#include <math.h>
#include <QTime>
#include <QDebug>
//#include <gsl_sf_trig.h>
//#include <gmp.h>
//#include <gmpxx.h>

//#define DEBUG_TAU_MOLE

void AtmosLayer::layerInit()
{
//	qDebug() << m_layerName << "layerInit()";

	m_lightInUp = new double[N_PI/2][2*N_PI][3];
	m_lightInD = new double[N_PI/2][2*N_PI][3];
	m_lightScatUp = new double[N_PI/2][2*N_PI][3];
	m_lightScatD = new double[N_PI/2][2*N_PI][3];
	m_lightScatUpTotal = new double[N_PI/2][2*N_PI][3];
	m_lightScatDTotal = new double[N_PI/2][2*N_PI][3];


	//initial out-scatter and incident light values
	for(int i = 0; i < N_PI/2; i++){
		for(int j = 0; j < 2*N_PI; j++){
			m_lightScatUp[i][j][0] = 0;
			m_lightScatUp[i][j][1] = 0;
			m_lightScatUp[i][j][2] = 0;

			m_lightScatD[i][j][0] = 0;
			m_lightScatD[i][j][1] = 0;
			m_lightScatD[i][j][2] = 0;

			m_lightInUp[i][j][0] = 0;
			m_lightInUp[i][j][1] = 0;
			m_lightInUp[i][j][2] = 0;

			m_lightInD[i][j][0] = 0;
			m_lightInD[i][j][1] = 0;
			m_lightInD[i][j][2] = 0;

			m_lightScatDTotal[i][j][0] = 0;
			m_lightScatDTotal[i][j][1] = 0;
			m_lightScatDTotal[i][j][2] = 0;

			m_lightScatUpTotal[i][j][0] = 0;
			m_lightScatUpTotal[i][j][1] = 0;
			m_lightScatUpTotal[i][j][2] = 0;
		}
	}

	if(m_zenithUpList.isEmpty()){
	    zenithAngleInit();
	    solidAngleCaculate();
	    scatAngleInit();
	}

	moleTauCaculate();
	m_tau = m_moleTau;
}


//caculate solid angle at NZ zenith angle
//m_solidAngle[i] is at Pi*i/N_PI,
//****Symetry: Solid_Angle(Pi*i/N_PI) == Solid_Angle(Pi-Pi*i/N_PI)***
void AtmosLayer::solidAngleCaculate()
{
//m_solidAngle[0] is solidAngle: zenith[0, ANGLE_STEP], azimuth ANGLE_STEP
//m_solidAngle[N_PI-1] is solidAngle: zenith[Pi/2-ANGLE_STEP, Pi/2], azimuth ANGLE_STEP
    for(int i = 0; i < N_PI/2; i++){
        m_solidAngle[i] = (cos(i*M_PI/N_PI) -
                   cos((i+1)*M_PI/N_PI)) * ANGLE_STEP;
//	Sin(iThetaI*M_PI/N_PI + ANGLE_STEP/2)*pow(M_PI/N_PI,2);

//qDebug() << i << m_solidAngle[i]; //tested 2011-10-31, night
	}
}

//#define SCAT_ANGLE_DEBUG
void AtmosLayer::scatAngleInit()
{
#ifdef SCAT_ANGLE_DEBUG
    qDebug() << "AtmosLayer::scatAngleInit()";
    QList<double> scatAngleList;
#endif
    m_scatAnglesSameDir = new double[N_ZENITH][N_ZENITH][2*N_PI];
    m_scatAnglesDiffDir = new double[N_ZENITH][N_ZENITH][2*N_PI];

    double muScat = 0;
    double thetaI, phiI, thetaS, phiS;
    thetaI = phiI = thetaS = phiS = 0;

    //same zenith direction
    for(int iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){
	thetaI = m_zenithUpList[iThetaI];
	for(int iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
	    thetaS = m_zenithUpList[iThetaS];
	    for(int iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){
		phiS = m_aziList[iPhiS];

		muScat = Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS);
		m_scatAnglesSameDir[iThetaI][iThetaS][iPhiS] = acos(muScat);
#ifdef SCAT_ANGLE_DEBUG
		if(!scatAngleList.contains(m_scatAnglesSameDir[iThetaI][iThetaS][iPhiS]))
		    scatAngleList << m_scatAnglesSameDir[iThetaI][iThetaS][iPhiS];
#endif
//		qDebug() << thetaI << thetaS << phiS << m_scatAnglesSameDir[iThetaI][iThetaS][iPhiS];
	    }
	}
    }
//    qDebug() << "AtmosLayer::scatAngleInit(), same dir over.";

    //different zenith direction
    for(int iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){
	thetaI = m_zenithUpList[iThetaI];
	for(int iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
	    thetaS = m_zenithDList[iThetaS]; //down, different
	    for(int iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){
		phiS = m_aziList[iPhiS];

		muScat = Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS);
		m_scatAnglesDiffDir[iThetaI][iThetaS][iPhiS] = acos(muScat);
#ifdef SCAT_ANGLE_DEBUG
		if(!scatAngleList.contains(m_scatAnglesDiffDir[iThetaI][iThetaS][iPhiS]))
		scatAngleList << m_scatAnglesDiffDir[iThetaI][iThetaS][iPhiS];
#endif
//		qDebug() << thetaI << thetaS << phiS << m_scatAnglesDiffDir[iThetaI][iThetaS][iPhiS];
	    }
	}
    }
#ifdef SCAT_ANGLE_DEBUG
    qDebug() << scatAngleList.size() << scatAngleList;
#endif
}

void AtmosLayer::zenithAngleInit()
{
    for(int iTheta = 0; iTheta < N_ZENITH; iTheta++){
	m_zenithUpList << (iTheta*M_PI/N_PI + ANGLE_STEP/2);
	m_zenithDList << (M_PI - (iTheta*M_PI/N_PI + ANGLE_STEP/2));
//		qDebug() << "zenith:" << iTheta << m_zenithUpList[iTheta] << m_zenithDList[iTheta];
    }
    for(int iPhi = 0; iPhi < 2*N_PI; iPhi++){
	m_aziList << (iPhi*M_PI/N_PI);
//		qDebug() << "azimuth:" << iPhi << m_aziList[iPhi];
    }
}

void AtmosLayer::setWaveLength(double wl)
{
	m_wl = wl;
}

QString AtmosLayer::getLayerName()
{
	return m_layerName;
}

void AtmosLayer::setTau(double tau)
{
    m_tau = tau;
    qDebug() << m_layerName << "set tau to" << m_tau;
}

void AtmosLayer::setResultSaveMode(char mode)
{
    m_saveMode = mode;
}

void AtmosLayer::setSavePath(QString path)
{
    m_resultPath = path;
    qDebug() << "Result gonner save at" << m_resultPath;
}

QString AtmosLayer::getResultPath()
{
    return m_resultPath;
}

//#define DEBUG_TAU_MOLE
void AtmosLayer::moleTauCaculate()
{
    if(m_hBottom == m_hTop){
	m_moleTau = 0;
	m_tau = 0;
	qDebug() << m_layerName << m_hBottom << m_hTop << "Molecule's Rayleigh scattering optical depth is:" << m_moleTau;
	return;
    }
    //initialize m_moleNumber, m_moleDensity
    if(m_moleNumber.isEmpty()){
	QFile moleFile("values/molecule_number_1976.txt");
	moleFile.open(QIODevice::ReadOnly);
	QByteArray moleFileData(moleFile.readAll());
	moleFile.close();

	QTextStream strm(moleFileData);
	while(!strm.atEnd()){
	    QString line(strm.readLine());
	    if(line.contains("#"))
		continue;
	    if(line.isEmpty())
		continue;
	    QTextStream strmLine(&line, QIODevice::ReadOnly);
	    QString altitude, numFromGround, density, numLayer;
	    strmLine >> altitude >> numFromGround >> density >> numLayer;

	    m_moleNumber.insert(altitude.toLong(),
				numFromGround.toDouble());
	    m_moleDensity.insert(altitude.toLong(),
				 density.toDouble());
	}
	//		qDebug() << m_moleDensity;
    }

/*
 --------   top molecule layer top, density(densTT)         ^  Z
							    ^
							    |
							    |
@@@@@@@@@@@  user input hTop                                |
							    |
							    |
							    |
 --------   top molecule layer bottom, density(densTB)      |
     *                                                      |
     *                                                      |
     *                                                      |
     *                                                      |
     *                                                      |
     *                                                      |
     *                                                      |
 --------   bottom molecule layer top, density(densBT)      |
							    |
@@@@@@@@@@@ user input hBottom                              | 0


 --------   bottom molecule layer bottom, density(densBB)
  */
    //find m_hBottom m_hTop corresponded layer in molecule data file
    QMap<long int, double>::const_iterator numB;
    QMap<long int, double>::const_iterator densityB;
    QMap<long int, double>::const_iterator numT;
    QMap<long int, double>::const_iterator densityT;

    numB = m_moleNumber.lowerBound(m_hBottom);
    densityB = m_moleDensity.lowerBound(m_hBottom);

    numT = (m_hTop > 1000000) ? (m_moleNumber.constEnd() - 1) : m_moleNumber.lowerBound(m_hTop);
    densityT = (m_hTop > 1000000) ? (m_moleDensity.constEnd() - 1) : (m_moleDensity.lowerBound(m_hTop));

    long hBB, hBT, hTB, hTT; //h-height, B-Bottom, T-Top
    double numBT, numTB;//num-molecule_number
    double densBB, densBT, densTB, densTT; //dens-density

    hBT = numB.key();
    numBT = numB.value();
    densBT = densityB.value();

    hTT = numT.key();
    densTT = densityT.value();

    numB--; densityB--; numT--, densityT--;

    hBB = numB.key();
    densBB = densityB.value();

    hTB = numT.key();
    numTB = numT.value();
    densTB = densityT.value();

    double moleNumLayer = 0;
    moleNumLayer += numTB - numBT;

    if(hBT != m_hBottom){ //when m_hBotton is right at boundery of molecule-file's layer, do not need extra add
	double densB = densBB + (densBT-densBB)*(m_hBottom-hBB) / (hBT-hBB);
	moleNumLayer += (densBT+densB)*(hBT-m_hBottom)/2;
    }
#ifdef DEBUG_TAU_MOLE
    qDebug() << "densB densBB densBT" << densB << densBB << densBT;
#endif
    double densT = densTB + (densTT-densTB)*(m_hTop-hTB) / (hTT-hTB);
    moleNumLayer += (densTB+densT)*(m_hTop-hTB)/2;

    m_moleTau = getMoleTauWholeAtmos() * moleNumLayer / m_moleNumber.value(1000000);
    qDebug() << m_layerName << m_hBottom << m_hTop << "Molecule's Rayleigh scattering optical depth is:" << m_moleTau;
#ifdef DEBUG_TAU_MOLE
    //tested 2011-11-16 night, bug fix for bottom.
    /*	AtmosLayer lay0("wholeAtmos", 0, 1000001, 0, 0);
 AtmosLayer lay10("wholeAtmos", 10, 1000000, 0, 0);
 AtmosLayer lay100("wholeAtmos", 100, 86000, 0, 0);
 AtmosLayer lay250("wholeAtmos", 250, 86000, 0, 0); */

    qDebug() << "wavelength" << m_wl <<
		"cross section" << sigma <<
		"m_moleTau" << m_moleTau <<
		"\tmoleNumLayer" << moleNumLayer <<
		"\nhBB" << hBB << "\tm_hBottom" << m_hBottom <<
		"\thTT" << hTT << "\tm_hTop" << m_hTop <<
		"moleNumBottomLayer" <<(densBT+densB)*(hBT-m_hBottom)/2<<
		"numBT" << numBT << "\n";

#endif

}

double AtmosLayer::getMoleTauWholeAtmos()
{
    QFile moleTauFile("values/tau_mole_Rayleigh.txt");
    if(!moleTauFile.open(QIODevice::ReadOnly)){
	qDebug() << "Error: cannot open file: values/tau_mole_Rayleigh.txt.";
    }
    QByteArray moleFileData(moleTauFile.readAll());
    moleTauFile.close();

    QTextStream strm(moleFileData);
    double moleTau = 0;
    while(!strm.atEnd()){
	QString line(strm.readLine());
	if(line.contains("#"))
	    continue;
	if(line.isEmpty())
	    continue;
	QTextStream strmLine(&line, QIODevice::ReadOnly);
	QString wl, tau;
	strmLine >> wl >> tau;

	if(m_wl == wl.toDouble()){
//	    qDebug() << "AtmosLayer::getMoleTauWholeAtmos(): wavelenght and tau: " << m_wl << tau.toDouble();
	    moleTau = tau.toDouble();
	}
    }
    return moleTau;
}

#define Log log
#include <math.h>
//tested ok 4-9, 2012, 22:56
QList<double> AtmosLayer::expHeightProfileCaculate(double ZStd,
						   double Zlower,
						   double Zupper,
						   int nlayers)
{
//    qDebug() << ZStd << Zlower << Zupper << nlayers;
    QList<double> heights;
    double z1 = Zlower;
    double z2 = Zlower;
    heights.append(z1);
    for(int i = 0; i < nlayers; i++){

	z2 = -1*(ZStd*Log((-1*exp(-1*Zlower/ZStd) + exp(-1*Zupper/ZStd))*(-1*(1/(exp(z1/ZStd)*(exp(-1*Zlower/ZStd) - exp(-1*Zupper/ZStd)))) + 1.0/nlayers)));
	heights.append(z2);
//	qDebug() << "expHeightProfileCaculate()" << i << z2;
	z1 = z2;
    }

//    qDebug() << "expHeightProfileCaculate(), heights:" << heights;

    return heights;
}
